Direct Numerical Demonstration of Sign-preserving Quasiparticle Interference via 
Impurity inside Vortex Core in Unconventional Superconductors 
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We perform large-scale numerical calculations self-consistently solving the Bogoliubov-de Gennes 
(BdG) equations in the magnetic field together with random impurities to directly demonstrate the 
typical quasi-particle interference (QPI) in the presence of vortices as observed by scanning tunneling 
microscopy /spectroscopy experiments in unconventional superconductors. The calculations reveal 
that vortex itself never works as a scatter causing the QPI pattern but vortex core containing 
impurity brings about the enhancement of the sign-preserving QPI peaks. Its origin is Andreev 
bound-states distorted by impurity, and all the measurement findings are consistently explained by 
the scenario based on the numerical results. 

PACS numbers: 74.20.Rp, 74.25.Op, 73.22.-f 
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I. INTRODUCTION 

Since the discovery of High-Tc cuprate superconduc- 
tors, the pairing mechanism has been a central issue in 
condensed matter physics over two decades. The su- 
perconducting gap symmetry was one of the most im- 
portant clues to identify the mechanism. Historically, 
the symmetry in High-Tc cuprate superconductors was 
proven to be d-wave by a direct observation of the half- 
fluxon in the tri-crystal junction and supported by other 
measurements^. Recently, a new type of symmetry, i.e., 
the sign-reversing s-wave pairing symmetry has been pro- 
posed as a candidate in brand-new iron-based High-Tc 
superconductors . 

The measurement of quasi-particlc interference 
(QPI) patterns by the scanning tunneling mi- 
croscopy/spectroscopy (STM/STS) is now listed as a 
powerful tool directly probing the pairing symmetry 3 - - — . 
The QPI pattern is obtained by Fourier-transforming 
the ratio Z{r,V) = g(r,V)/g(r,—V) of the conductance 
maps g(r, V) = dI/dV(r 7 V) for the position r and the 
bias voltage ±V. So far, it has been widely accepted that 
a relative sign difference inside superconducting order 
parameter can be directly detected by the magnetic field 
dependence of its Fourier transformation \Z(q,u>)\. The 
reason is that if a vortex works as a magnetic impurity 
then the number of magnetic impurities increases with 
increasing the magnetic field and the sign-preserving 
scattering 5 dominates over the sign-reversing one. How- 
ever, it is just an intuitive idea, and there is no direct 
theoretical confirmation except for an indirect support 
through a perturbative approach?. Generally, vortex 
breaks the translational symmetry and makes low-energy 
Andreev bound states inside the core. Thus, the vortex 
is never a simple impurity and beyond the reach of any 
perturbative approaches. In iron-based superconductors, 
the QPI experimental result^ is presently regarded as 
one of a few crucial facts supporting the sign-reversing 
s-wave symmetry. We have to rush to confirm the 
intuitive assumption. 



In this paper, we therefore investigate QPI in the pres- 
ence of vortices through a direct numerical calculation on 
the self-consistent BdG equations. A main finding of the 
present paper is that vortex works as a scatter break- 
ing the time-reversal symmetry only when impurity is 
captured by the core. Thus, this finding indicates that 
vortex pinning, i.e., so-called the core pinning is crucial 
in supporting the previous QPI scenario in the magnetic 
field. 

To numerically obtain the QPI map, we need to per- 
form a large-scale calculation in real-space, since the QPI 
map is obtained by Fourier transforming a sufficiently 
wide real-space data. However, the computer resource to 
numerically solve the BdG equations becomes too huge, 
as long as one employs the conventional Hamiltonian- 
matrix diagonalization-scheme whose computational cost 
increases in A~ 3 manner, where N is the number of grid 
points in real-space. In this paper, we instead imple- 
ment the Chebyshev-polynomial expansion schemed— 
as a self-consistent solver of the BdG equations^. We 
stress that the present real-space area is wider than 10- 
nm scale square, which is much beyond the conventional 
size. 



II. MODEL AND METHOD 



The target tight-binding model Hamiltonian H 



%BCS + %e 



Here, BCS Hamiltonian, Hbcs 



h.c. 



where 



- T,ij,a(tij + ^ii) c i* c 3° + Ei 

c\ a creates an electron with spin a at site i and jj, denotes 
the chemical potential. The hopping integrals tij include 

the Peierls phase factor ty = exp J^ J A(r) ■ dr 

in the presence of the magnetic field, where A(r) is the 
vector potential and <fio = hc/2e is the flux quantum. 
The impurity part of the Hamiltonian can be written as 



H; 



imp 



E, 



y imp c f c- 



One diagonalizes 1-L to solve 
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the BdG equations written as 



K A- ■ 

At 



-K 



>-j 



u a (rj) 
v a {rj) 



— E a 



u a (ri) 
v a (ri) 



(1) 



Here, K id = - ( j u-'V? inp )<5j 7 - , and Ajj = A(r i ,r i ) = 

Vij X^q^ u a{ r j) v a ( r i)f(Eu), where N is the number of 
the lattice sites, denotes the pairing interaction, and 
f(x) is the Fermi distribution function. In this pa- 
per, the hopping is restricted only in the nearest neigh- 
bor (tij = t) for simplicity. The chemical potential 
fj, = — 1.5t, and the pairing interaction is given only 
on the link of the nearest neighbor, Vij = V = —2.2t. 
We self-consistently calculate d x 2_ y 2-wave order param- 
eter, A d (n) = ViA^i + A^ ti - Ay, t - A_ 5i i)/4 with 



A ±i>i = A(r,,n ± e)exp [i^ £ i+ri±S)/2 A(r) ■ dr 
where x and y denote the unit vectors in a square lattice. 

Let us briefly show how to solve self-consistently 
the BdG equations JT]) including the gap equa- 
tion with use of the Chebyshev-polynomial expansion 
scheme. The mean-field can be expressed as (ciCj) — 

-^JZo du f(u)e(j) T d{u)h[i), wner e [e(i)] 7 = S hl , 
and [/i(i)]-y = <5i+Ar l7 . The spectral density d(u) is given 
by G R (u;) — G a (oj) whose Dirac's delta functions are ex- 
panded by a series of Chebyshev-polynomials 1 -. Then, 
one can rewrite the gap equation at zero-tcmpcrature as 
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A(n, rj ) 



2Vn 



sin[n arccos(— b/a)] 



n=l 



(2) 



where n is the order of the Chebyshev-polynomial, and n c 
denotes a cutoff parameter. h n (i) is calculated by the re- 
currence formula, h n+ i{i) = 2(H — lb)/ah n (i) — h n -\(i) 
and ho(i) — h(i) and h\(i) = 2(71 — lb)/ah(i), where 
the renormalized factors a and b are set in the order of 
the band- width, e.g., b = — fi, respectively. We point 
out that calculation results are insensitive to choice of 
these parameter o 17 ! 18 . This drastically reduces the self- 
consistent calculation costs since the calculation algo- 
rithm is perfectly free from the use of any diagonalization 
scheme. We always set n c = 1000 through the present 
calculations, which is confirmed to be enough 1 - 8 . Since 
each grid point is completely independent in the vector 
formula, we can efficiently solve these equations by using 
parallel multi-core computers. 

The calculation target is 160 x 160 square lattice, which 
conventionally requires a full diagonalization of 51200 x 
51200 matrix. If the lattice constant ~ A, then its real- 
scale area is beyond 10 nm square. In this case, it is 
a quite hard task for the conventional diagonalization 
scheme to self-consistently solve the problem, while the 
present scheme can finish the self-consistent calculation 
during about 10 minutes when using 128 cores. 
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FIG. 1. (Color online) (a): A spatial profile of the d x 2_ y 2- 
wave order-parameter amplitude |A,j(rs)|. (b): A q-space 
profile of Z(q,E), which is obtained by Fourier transforming 
the quantity Z{r,uf) = N(r,cj)/N(r,—u}) at the energy cj = 
0.04i. The system size is 160 x 160 square grids with 10 
randomly distributed impurities without vortex. 



III. RESULTS 

Here, let us show how to demonstrate that QPI in the 
presence of the magnetic field actually detects the sign- 
preserving quasi-particle scattering. First, we examine 
the quasi-particle scattering by non-magnetic impurities 
without vortices. Second, we study QPI for a vortex 
lattice without any impurities to know whether vortex 
works as a magnetic impurity. Third and fourth, we in- 
vestigate two vortex systems with random impurities. In 
the third case, all impurities are away from any vortex 
cores. On the other hand, one impurity locates inside the 
vortex core in the forth case. Through the comparison 
between the third and fourth cases, we find that impu- 
rity locating inside the vortex core is essential to identify 
the relative phase difference of the superconducting order 
parameter. 

At first, let us examine the QPI with no vortex. We 
introduce randomly-distributed impurities as V^ lmp = 
0.3tS(r - r™? 4 ) (i = 1, • • • 10). As shown in Fig. H^a), 
the order parameter amplitude |Ad(rj)| shows relatively 
small values only around the impurity sites. In this case, 
as seen in Figfljb), the quantity Z(q, E), i.e., the Fourier 
transformation of Z(r,u)(= N (r , uj) / TV (r , —to)) is equiv- 
alent to those in the previous theoretical studies using 
the diagrammatic approach^. This map is regarded as a 
typical QPI pattern, which clearly reflects the scattering 
feature of quasi-particles on the Fermi surface (see the 
left panel of Fig. 4 for the present Fermi surface) . 

Second, we study a vortex lattice without any im- 
purities. According to the standard way introducing 
vortice a 14 ' 19 , we make a square vortex lattice, whose unit 
vector is in the direction of 45° from a-axis of the original 
tight-binding model (see the inset of Fig. [5]). Then, the 
local density of states (LDOS), N(r, E) is calculated and 
the gap edge is found to be ~ 0.252; as seen in the left 
panel of Fig. [5] The QPI map is displayed in the right 
panel of Fig. [51 where we note that the map is irrelevant 
to QPI. Namely, one can not find out any characteristic 
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FIG. 2. (Color online) Left panel: Local densities of states 
far from a vortex core and inside a vortex core without any 
impurities. Inset: A spatial profile of the d a .2_ y 2-wave order- 
parameter amplitude |Ad(ri)|. Right panel: A q-space profile 
of Z(q, E). For the definition Z(q, E) and the system size, see 
the caption of Fig. Q] 



wave-vector peaks arising from the impurity scattering 
in contrast to Fig. [TJ In the case, Andreev bound states 
inside the vortex core just produce such a pattern only 
around q ~ 0. Thus, it is found that vortex can not be 
regarded as a magnetic impurity scatter. 

Now, let us introduce random impurities into the vor- 
tex lattice. As the third case, 10 impurities are ran- 
domly distributed, but all impurities are away from any 
vortex cores. In this case, their impurity potential is 
not so deep and their density is so dilute that the vor- 
tex lattice is not almost distorted compared to the case 
without impurities. Such a situation is called "floating 
lattice" , in which an energy gain by the vortex lattice for- 
mation dominates over that by the vortex core pinning. 
As shown in Fig. EJb), the QPI pattern in the floating 
lattice is almost regarded as a sum of QPI's in the first 
and second cases. 

Here, let us move each impurity site into the vortex 
core in order to create the so-called " pinned lattice" . In 
this case, we find that the calculated Free energy is re- 
duced compared to that of floating lattice. As the fourth 
case, we move an impurity marked by the red line cir- 
cle in Fig. [3ja) to the marked vortex core as shown in 
Fig. [3Jc). Then, the marked vortex of Fig. (3Jc) is as- 
signed to "pinned vortex", and the QPI in the presence 
of the pinned vortex is displayed in Fig. [3Jd). By com- 
paring Fig. [3^d) with Fig. [3jb), one finds that intensified 
points are located close to yi = (ir, 0), q% ~ (2.3, 2.3) and 
qs ~ (1.2, 7r). As these q's are displayed in the left panel 
of Fig. 21 qi, q 2 and qs, correspond to the scattering vec- 
tors represented by k F (0) - k' F (n), fc£(7r/4) — fc^(57r/4), 

and k F (3w/8) — A;^(ll7r/8), respectively, where k F (8) 
and k F (9) are the Fermi wave vectors after and before the 
scattering, and denotes the position angle on the Fermi 
surface. One can finds that qi and 93 vectors mean sign- 
preserving scatterings. On the other hand, one notices 
for g2-vector that one needs a numerical check in judg- 
ing whether qi really corresponds to a sign-preserving 
scattering or not. The result is displayed in the right 
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FIG. 3. (Color online) (a) and (c): Spatial profiles of the 
d x 2_ y 2-w&ve order-parameter amplitude |Ad(ri)|. (b) and 
(d): q-space profiles of Z(q,E). For the definition Z(q,E), 
see the caption of Fig. [I] We consider the 160 x 160 square 
lattice system with 10 randomly distributed impurities with 
vortices, (a) and (b): Any impurities are not located inside 
vortices ("unpinned vortex"), (c) and (d): One impurity is 
located inside a vortex ("pinned vortex'). 



panel of Fig. |H which distinguishes the sign-preserving 
with the reversing one in q-space quadrant by actually 
examining q- vectors around the Fermi surface (±0.1i). 
From the right panel of Fig. [4] it is found that q-i is also a 
sign-preserving vector. These facts indicate that the sign- 
preserving scatterings are intensified by the emergence of 
the pinned vortex core. In addition, we should note that 
the vortex arrangement influenced by the vortex pinning 
is not important for QPI phenomena, since these origins 
are localized near a vortex-core and an impurity. We also 
note that these QPI phenomena are not originated from 
the disordered effects since we put only ten impurities 
whose density is 0.04% in our system, which can be as- 
signed as a dilute-impurity doped case. In reality, the 
number of pinned vortex cores increases when increasing 
the magnetic field. Then, one can expect that only the 
sign-preserving q-points grow in the QPI pattern. We 
emphasize that it is consistent with the measurement re- 
sults. The presence of pinned vortices is essential for 
the phase sensitive detection of the superconducting or- 
der parameter. Indeed, such sign-preserving scattering 
peaks are non-observable in too clean sample 20 . 



IV. DISCUSSION 

Finally, let us discuss the present numerical scheme 
and the obtained QPI results. The first issue is an advan- 
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FIG. 4. (Color online) Left panel: A fe-space profile of the 
Fermi surface (red solid curve) in the square lattice model 
with fj, = — 1.5t and the white and shaded areas represent 
signs of d-wave superconducting gap on fc-space. The arrows 
denote typical sign-preserving scattering q- vectors (See the 
text). Right panel: A q-space profile of the sign preserving 
(+) and reversing (-) on the superconducting gap. 



tage of the Chebyshev-polynomial expansion scheme in 
self-consistently calculating the BdG equations. In order 
to obtain a QPI map by Fourier transforming a real-space 
conductance map, one needs to solve the BdG equations 
on a sufficiently large real-space. Thus, there has been 
so far no study to directly calculate QPI's. This clearly 
indicates that the conventional diagonalization scheme is 
never a practical way in studies like the present topic 
and an alternative more efficient scheme is demanded. 
Very recently, Covaci et al, have proposed a quite ef- 
ficient scheme based on the Chebyshev-polynomial ex- 
pansion to examine inhomogeneous superconductivity^. 
The main idea is an expansion of the Green's function by 
a set of Chebyshev-polynomials, which guarantees a dras- 
tic reduction of the calculation cost. However, Covaci et 
al., left a self-consistent calculation as the future work. 
Thus, the present work is the first self-consistent large- 
scale BdG calculation including the vector potential. We 
confirm in the present case that the computational cost 
totally scales with 0(N 2 ) on the real-space grid size N 
and find that the number of the self-consistent iterations 
is almost equivalent with that in the matrix diagonaliza- 
tion. This means that the calculation cost is reduced to 
~ 1/N compared to the diagonalization scheme. More- 
over, we developed several techniques in order to acceler- 
ate the calculation (see RefJ£ for the details). The sec- 
ond issue is an origin of the present QPI pattern in the 
presence of the pinned vortex. Previously, the resonant 
Andreev scattering has been regarded as the scattering 



mechanism around a vortej*2ii£. However, our numerical 
calculations reveal that the spatial variation of the order 
parameter induced by a vortex does not work as an ef- 
fective scatter. In fact, as seen in Fig. [51 we can not find 
out any characteristic wave- vector peaks unlike the typi- 
cal QPI map as Fig.[TJ Thus, the QPI origin of the pinned 
lattice is explained as follows. The Andreev bound states 
around a vortex core are slightly distorted by impurity 
inside the vortex core^i. A set of such distorted states 
scatter quasi-particles like an impurity while still keeps 
the angular momentum breaking the time reversal sym- 
metry. Thus, we easily notice that a pinned vortex be- 
comes a sign-preserving scatter and the number of such 
scatters increases with increasing the magnetic field. On 
the other hand, the number of the sign non-preserving 
scatter contrarily decreases with increasing the magnetic 
field, since the number of the scalar impurities outside the 
vortex core decreases with increasing the magnetic field. 
These consequences are consistent with the experimental 
results. In addition, our discussion can be easily applica- 
ble to the iron-based superconductors^, since the present 
above scenario is considered to be also kept in multi-band 
superconductors. The realistic calculation for iron-based 
superconductors is a future work. 



V. CONCLUSION 

In conclusion, we self-consistently performed large- 
scale BdG calculations in the presence of vortices to- 
gether with random impurities by using the Chebyshev- 
polynomial expansion scheme to investigate QPI in the 
magnetic field. Our calculations concluded that the a 
vortex core distorted by impurity works as an impu- 
rity breaking the time reversal symmetry. The micro- 
scopic finding well explains the observed experimental 
data, e.g., the QPI peaks relevant to the sign-preserving 
quasi-particle scattering grows while the sign-reversing 
one diminishes with increasing the magnetic field. Our 
direct calculations confirmed that the magnetic field de- 
pendence of QPI is a true phase-sensitive tool for uncon- 
ventional superconductors. 

We would like to thank Y. Ota and R. Igarashi for help- 
ful discussions. The calculations have been performed 
using the supercomputing system PRIMERGY BX900 
in Japan Atomic Energy Agency. 
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